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ABSTRACT 

In  section  A,   the  theoretical  basis    for  representing  a  rational 
number  losing  the  prime    factorization   of  its    minimal  integer   ratio   is    given, 
and  the  resulting  algorithms    for  performing  addition,   subtraction,   multi- 
plication,  and  division   are   developed.      The   conjecture   that  these    algorithms 
will  normally  require   fewer  machine   operations   than   algorithms  based  on 
storing  rational  numbers    as   two   long  integers   is    supported  by   a  worst   case 
machine   operation   count  estimate   for  each.      These   operation   counts   are 
applied  to  parameters    derived  from  test   results. 

Section  B  is    a  self-contained  user's   manual   for  RATION,    a 
rational   arithmetic  package  written   in   FORTRAN    and  based  on  the   above 
algorithms . 


Digitized  by  the  Internet  Archive 
in  2013 


http://archive.org/details/rationarationala572brow 


11 
TABLE   OF   CONTENTS 

Page 

A.  RATIONAL  ARITHMETIC  ALGORITHMS 1 

1.  Theoretical  Foundations 1 

2.  Conversion  Algorithms 1+ 

3.  Arithmetic  Algorithms    7 

k.      Package  Organization 11 

B.  RATION— USER'S   MANUAL 17 

5 .  User  Program 17 

6.  I/O   and  Utility  Instructions 18 

7.  Arithmetic   Instructions    20 

LIST  OF  REFERENCES 22 

APPENDICES 

I.      Sample  Program 23 

II.      RATION   Listing 25 


A.      RATIONAL  ARITHMETIC  ALGORITHMS 

1.      Theoretical  Foundations 

A  package    of  FORTRAN   subroutines   has  been  programmed  to   perform 
arithmetic   operations   on   a  special  representation   for  rational  numbers 
using  algorithms   that   attempt   to  minimize    computer  execution  time   at   the 
expense   of  space   required  for  the   storage   of  the  operands    and  the  program 
itself.      The   algorithms    and  representations    depend  on   some  elementary 
properties   of  the  sets   of  integers    and  rational  numbers    as    developed  in 
number  theory.      These   properties    are    defined  below.      See    [l]    for    fuller 
explanations . 

Definition   1:      The  set  I   of  integers    consists   of  the   counting   (natural) 
numbers    (1,2,3,...)    and  their  negatives    (additive   inverses) 
and  0 . 

Definition  2:      R  is    an  element   of  the  set  R  of  rational  numbers   if  it    is 
any   rati  or/ 1    /I  _,   I      e   I_,   I  _   e   I_.      Note   that  this    defines    an 
equivalence    class    for  each   rational  number.      The   rational  number 
is   the  equivalence    class    (I    *I    )/(l    *I    )    for  I      e    I.      If 
Ip   >    0   and  there   exists  no   integer  I,    such   that   I*    divides  both 
I      and  I    ,   then   I   /I      is   the  minimal   element  of  the  equivalence 
class    and  will  be   referred  to   as   the  rational  number  R  =   I    /I    . 
See    [6]    for  further  explanation. 

Definition   3:      A  prime  number  is    any  positive   integer  which  is    divisible 

without   remainder  by  exactly  two   distinct  integers — 1   and  itself. 

Definition    k:      Any   integer  I    can  be   uniquely   represented  by  the  product   of 
some  set  of  prime  numbers.      This    unique   set  is    called  the  prime 
factorization   of  I. 


Using  the   above    concepts  ,   one    can   see  that   any  rational  number 
involves   two   integers — each  of  which  has    a  unique  prime    factorization. 
The  prime    factorization  has    a  unique  representation  when  placed  in  the   form 

h  =  Pl     P2    —Pn 

a.    >   0 

1 

where  p     is   the  i-th  ordered  distinct  prime   appearing  in   the   factorization 
1 

and  a.    is   the  number  of  times   it   appears. 

i 

Since   a.    >    0,    a  representation   for  l/ip   is 

~"61     ^2        T3m 
P-L       P2      ...Pm 

e.  >  o. 

i 

Since   the  rational  number  R  is   I    /I    ,    a  representation   for  R  is 

a1-g1    a2-§2        ft  -g 

Pl  P  --'Pq 

where    {p.}   is   a  reordering  to  include   all  elements   of  {p.},   (p.),    and 

ft.    =   0   if  p.    A    {p. } 
1  1   T        1 

3.    =   0   if  p.    £    {p.}. 
i  i   T        l 

This   representation   can  be   used  in   a  computer  by  storing  a  string 

of  nodes,   each  node    containing  the  exponent   of   a  prime   associated  with  the 

node.      If  we   let   an   indexed  array  P(l)    contain    our  primes  with  P(l)    equal 

to  the  total  number  of  primes    defined  in  the   array    (this   may  be   less  than 

its    dimension)    and  P(2)    =  2 ,   P(  3)   =    3,   P(U)    =   5 ,  P(  5)   =   7,. .  .  ,  then  each  node 

will  take   the   form  of   (#p.    )   where  y.    is  the  exponent   a.-0.    of  the  prime  p. 

with   index  #p.  . 
l 


One   additional  node   is  needed  to  specify  the  sign   of  the  rational 
number  and  the  number  of  nodes   occurring  in  this  particular  rational 
number's   representation  since   this  will  seldom  be   the   full  length  of  the 
array.      By  using  a  FORTRAN   array  dimensioned  R(2,K)  ,   10   <_  K   <_  25,  we    can 
set   R(l,l)    so 

signum[R(l,l)  ]   =  signum[R] 
|R(l,l)|    =  number  of  nodes   including  the  bookkeeping  node. 
Also  set  R(2,l)   =   index  of  largest  prime   in   the   array,   R(2,J)   =   index  of 
prime  number  taken  to  the  R(l,J)   power   (J  >_  2).      In  keeping  with  our 
notation,   -1,   0,   1  have   R(l,l)    =   -1,   0,   1,    respectively,    and  R(2,l)    =   1 
by   convention. 

Henceforth,   a  rational  number  R(2,K)    in   this    form  will  be   called 
the  indexed  prime    factorization   representation    (factored  form  for  short) 
of  the  rational  number  R.      It  will  be  written   in  the   form  used  in   the 
following  examples . 

Ex.    1: 


Ex.    2 


Ex.    3: 

A  representation  will   also  be  needed  for  very  long  integers   if 
numerators    are  to  be   added  or  subtracted.      If  a  suitable  base   is    chosen  such 
that   any  number  multiplying  a  digit   of  the  integer  will  leave  the  product 

completely  within   one   computer  word ,   thus   retaining  the   carry,   an   integer   can 

*  r-1   if  A  <   0 

signum[A]   =  /    0   if  A  =   0 

I  1   if  A  >   0 


1* 


be   represented  as    a  string  of   computer  words,   each  word  representing  an 
(extremely  large)    digit.      Again,    a  bookkeeping  entry  is  needed  to  specify 
the  number   of  digits    defined  and  the  sign   of  the   integer  I.      Otherwise, 
all  entries   of  the   array  would  have  to  be    zeroed  before   a   change   is  made. 
Therefore,    |P(l)|    is   set  to  the  length  of  the   array,   signum[P(l)  ]   to 

signum[l],    and  P(j)   to   the   coefficients   of  BASE**(J-2)    for  J  >_  2. 

The   following  examples   illustrate   the   representation   of  long  integers. 

The  base   is   10000   so  the   digits   0000   to   9999   are   represent  able. 

Ex.    k:  U7926U3  =    (3,261+3,i+T9) 

Ex.    5:  -177133=   (-3,7133,17) 

-1  =    (-2,1) 
Ex.    6:  0  =    (2,0) 

1  =    (2,1) 

2.      Conversion  Algorithms 

It  will  be   shown   in   section   3  that   during  addition,   factored 
forms   of  two  integers   must  be   converted  to   long  integer   forms,    added,    and 
the  resulting  long  integer   converted  back  to   the   factored  form.      Two 
algorithms    are  necessary   for  these    conversions. 

Consider  the    case   of  a  factored  form  integer  A   (all  positive 
exponents)    and  a  list   of  prime  numbers  none    of  which  exceeds   one    computer 
word  in   length.      This   is   reasonable  since  n(x) — the  number  of  prime  numbers 

less   than   x — is    approximated  by  x/log     x   [2]    and  on   any  IBM  360    computer 

0 
this   allows   the  program  to   allocate   10     prime  numbers — much  more   than   could 

ever  be   stored  effectively.      Then   if  long  integer  II  is   set   to  1  and  each 

element  Il(2),...    is   multiplied  by  prime  P(J)    (carrying  properly)    and  each 

P(J)    is   multiplied  times   II   a.    times    as   J  is   stepped  through   all  primes  with 

J 

positive   a.    in  A,   II  is    converted  to   a  long  integer   form  of  A.      This   is 

J 

shown   in  Algorithm  C,  which,   like   all   algorithms  herein,   is   expressed  in 
Knuth's  notation    [3]. 


ALGORITHM  C(A,Il) :      Convert   factored  form 
integer  A  to   long  integer  II 


CI.  [Initialize.  ]   Il(l)  +■  2,   Il(2)  -*-  1,   M  *-   2. 

C2.  [Prime   factor  loop.]     J  •*■  2. 

C3.  [Exponent   loop.  ]     K  «-  1. 

Ck.  [Long  integer  multiply.  ]      L  *■  2,   Il(L+l)   =   0. 

C5.  [Prime  P  times   II.]      IL(L)  <-   IL(L)*P(A(2  ,J) ) . 

If  I1(L)    >   BASE,    I1(L)  «-   I1(L)    -    LI1(L)/BASEJ*BASE, 

Il(L+l)   «-    LH(L)/BASEJ. 

If  L  4  M,  then   L +-    L+l ,    go  to  C5. 

If  Il(L+l)   ^   0,   M^    M+l. 

C6.  [End  loop.  ]      If  K  4  A(l,j)  ,  K  *■   K+l ,    go   to   CU. 
CT.  [End  loop.  ]      It  J  4    |  A(l,l)|  ,  J  «-   J+l ,    go  to   C3. 
C8.  [Set  bookkeeping  node.]      Il(l)  <-   M*signum(A(l  ,l) ) . 
END. 

Going  the  other  way   (converting  a  long  integer  to   factored  form) 
one   must    check  if  each  prime   divides   the  integer,    and  if  so,   add  1  to  the 
exponent  of  that  prime   in  the   factored  form  result;   if  not,   try  the  next 
prime  stopping  when   II   =   1.      This   is   shown   in  Algorithm  G.      It  is   assumed 
that  a  function  IDIV(ll,P)    for  dividing  a  long  integer  II  "by   a  single 
word  integer  P  is   available.      This  would  be   similar  to  the  multiply   routine 
in   C3-C6  of  Algorithm  C. 


ALGORITHM  G(I2,A):      Convert   long  integer 
II  to   factored  form  A 


Gl.  [Initialize.  ]      M^2,I^2. 

G2.  [Check   for  II  =   1.  ]      If  Il(l)   =   2   and  Il(2)   =   1,   then 
A(l,l)  <r.   I*signum   (ll(l)),  A(2,l)  +"   M,   END. 

G3.  [Check  remainder  on   division  by   a  prime.] 

TEMP  =   IDIV(I1,P(M) ) .      If  there  is    a  remainder, 
M  <-   M+l,    go   to   G2;    otherwise,    II  ■«-   TEMP. 

Gk.  [Increment   factor   count.]      If  M  4  A(2,l),  then 
I  ■*■   1+1,   A(2,I)  <■  M,   A(1,I)  +■    0.      Set 
A(l,l)«-    A(1,I)   +   1.      Go   to   G2. 


Algorithm  G  must    access    all  the  primes    sequentially  so  the 
factored  form  has  the  primes   in   order.      Algorithm  C   accesses   the  primes 
sequentially  "but    could  "be  written  to   access  them  randomly.      If  the   list 
of  primes   is   needed  sequentially,   a  sequential   access    device    could  easily 
be  used  to  store  the   list  ,   although  the   current  program  package  stores 
them  in   core.      With  these    conversion   algorithms    defined  the   arithmetic 
algorithms    can  "be   completely  specified. 

The   method  for   choosing  the  base    for  long  integers    can  now  be 
defined  also.      Since   any  digit  may  be  multiplied  by  the   largest  prime   in 
the  list,   this  product   must   still   fit  into   one   signed  integer  on  the 
computer  in  use  to  allow   carries.      Thus,  BASE   <_  MAX/P(P(l)) — where  MAX  is 
the    largest   signed  integer  representable   on  the   computer  in   use — and 
P(P(l))    is  the  largest  prime  number  that  will  ever  occur  in  the    list   of 
primes  P(l) . 
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3.      Arithmetic  Algorithms 

Four  operations    are    closed  on  the  set   R_  of  rational  numbers: 
addition,   subtraction,   multiplication,    and  division    [k].      All  four  of 
these   can  be   implemented  easily   and  effectively  with  the  previously 
described  representations    for  rational  numbers    and  integers.      In   fact, 
multiplication  and  division  become  trivial  for   factored  form  rationals. 

To   show  this,    consider  two   factored  form  rationals 

Since   their  product   is   the  sum  of  exponents   of  common   factors  ,   it   can  be 

found  simply  be   setting  a  marker  IA,   IB    for  each    and   advancing  it    through 

aiA  3IB 

(#pTA),    (#p      ).      At   each   configuration   reached,   either  #pT.    >    #P-,-n   or 
1A  lis  IA  IJd 

y± 

#p        >   #pT,    or  #pT.    =   ^P-r-r,.      In  each    case,   the  node   in   the   result,    ( #p      ) 

Ixi  IA  IA  1x5  1L 

which  has   its   own  marker  IC ,   takes   on  the   lower  indexed  prime   and  the 
corresponding  exponent    and  increases  the   corresponding  marker  by  1.      It 
also   increases   its   own  marker  by  1 .      If  the  prime   index  on  both  operands 
is  the  same,   their  exponents    are    added.      In  us  ,   for  the   configuration   above 

YIC  =    3IB'  YIC   =    aiA'    °r  YIC   =    aiA  +   ^IB'   resPectively-      Tnis   leads   to 
Algorithm  M. 


ALGORITHM  M(A,B,C)  :      C  *-   A*B 


MO.  [Set  multiply   indicator.]      IN  *   1. 

Ml.  [Initialize.  ]      IA  *   2,    IB*-   2,    IC  *-   2, 
II  *   MAX(A(2,l)  ,   B(2,l))  ,   I  *-   2. 

M2.  [Search  for  next  larger  prime.]  If  A(2,IA)  =  B(2,IB)  =  I, 
go  to  M5.  If  A(2,IA)  =  I,  go  to  M3.  If  B(2,IB)  =  I ,  go 
to  Mh.      Otherwise,   I  *   1+1. 

M3.  [A  factors    C.  ]      C(2,IC)  *   I,    C(l,IC)  *■   A(l,IA), 
IA  *    IA+1,    IC  *   IC+1.      Go  to  M6. 

Mk.[B   factors    C.  ]      C(2,IC)  *    I,    C(l,IC)  *    B(l,IB)*IN, 
IB  *   IB+1,    IC  *   IC+1,    go   to  M6. 

M5.  [A  and  B   factor  C.  ]      C(2,IC)  *   I,    C(l,IC)  *   A(l,IA)   + 
B(l,IB)*lN.      IA  *   IA+1,   IB*    IB+1,   IC  *   IC+1. 

M6.  [Loop  end.]      If  I  4  II,  I*-   1+1,   go  to  M2;   otherwise, 
if  C(1,I)   =   0,1*   1-1,   go  to  M6. 

M7. [Bookkeeping  node.]      C(l,l)  *   I*signum(A( l,l)*B(l,l) ) . 
C(2,l)  *    C(2,I).      END. 


For   division,    C  =  A/B— the  same   method  of   advancing  through  the 
nodes   of  A  and  B   and  comparing  the  indices   of  the  primes   holds  ,   except 
6IB  is   rePlaced  by   -3IB  since   division   is   accomplished  by  subtracting 
exponents   of  common   factors.      This   leads  to  Algorithm  D,  which   sets   to  -1 
a  marker  to  be   multiplied  by  3      . 

ALGORITHM  D( A, B,C)  :      C  *-    A/B 

Dl.  [Set   divide   indicator.]      IN*   -1. 
D2.Go  to  ML. 


For   addition,   a  least    common   denominator   (l.c.d.  )    of  the  two 
rational  operands   must  be   found.      This   is    actually  the   least    common 
multiple    (l.c.m.  )    of  the  two   denominators.      The   denominators   of  the 
operands   can  he   found  easily  since  they  are   representable  by  extracting 
only  the  nodes  with  negative   exponents.      Their  l.c.m.    can  be   found  by 
choosing  the   more  negative   of  any  negative   exponents    for  any  primes   in 
either  of  the  two   factored  form  operands.      Taking  absolute  values    gives 
the   least   common   denominator.      This    result   is   shown  in  Algorithm  L. 


ALGORITHM  L(A,B,R):      R^    least    common 
denominator  of  A   and  B 


LI.  [Initialize.  ]      IA^2,IB^2,IR^2,I^2, 
II  +■   MAX(A(2,1),   B(2,l)). 

12.  [Search   for  next  prime.]      If  A(l,IA)    =   B(l,IB)   =   I, 
go  to  L5.      If  A(1,IA)   =   I,   go  to   L3. 
If  B(1,IB)   =   I,    go   to  LU. 

L3.  [A   factors    R.  ]      R(l,IR)   ■*■    -MIN(A(l,IA)  ,0)  . 

IA<-   IA+1.       If   R(1,IR)    4   0,   R(2,IR)«-    I ,   IR  ■*■   IR+1. 
Go  to   L6. 

LU.[B    factors   R.  ]      R(l,IR)  ■*■   -MIN  (B(l  ,IB)  ,0  ) . 

IB  *■   IB+1.      If  R(l  ,IR)    4   0  ,  R(2  ,IR)  «-    I ,  IR  «-    IR+1. 
Go  to   L6. 

L5.  [A  and  B    factor  R.  ]      R(l,IR)    =   -MIN(A(l,IA)  ,    B(l,IB)  ,0). 

IA  «-   IA+1,    IB  *    IB+1.      If  R(1,IR)    4   0,   R(2,IR)  +•   1 9   IR +■    IR+1. 

L6.  [Bookkeeping  node.]      If  I  ^  II  ,   I  ^   1+1 ,    go  to  L2  ; 
otherwise,   R(l,l)  +   IR,   R(2,l)  *■   R(2,IR). 
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If  D  represents  the  l.c.d.  ,   then   for  A,B   rational  numbers 

A*D,  B*D  are   factored  form  integers.      These   integers    can  be    converted 

to  long  integers  "by  Algorithm  C.      If  each  succeeding  element  of  the 

corresponding  integers    IA(K)  ,  K  =  2(l)|lA(l)|  ,   IB(L),    L  =   2(l)|lB(l)| 

is    added  with   appropriate    carries   or  borrows    as   in  Algorithm  I,   the 

numerator  can  be   found.      Note   that 

,       ,    v    _  /0   if  x  <   0 
step(x)   -{lifx70. 

Converting  the  numerator  to   factored  form  using  Algorithm  G  and  then 
dividing  by  D,  the    denominator,   using  Algorithm  D,    gives   the   result. 
Since  the  signs    of  the   operands    do  not  matter  once   they   are    attached 
to  the  long  form  integers  ,   subtraction   can  be    carried  out  by   changing  the 
sign   of  the  subtrahend  before   and  after   addition.      Thus,   addition   and 
subtraction   are  performed  as    in  Algorithms   A  and  S. 

ALGORITHM   1(11,12):       II  -e    11+12 


11.  [Initialize.  ]      IP1  «-    |  Il(l)  |  ,   IP2  «-    |l2(l)|, 

I  «-    2,    II  <-  MAX(IP1,IP2)  ,    CNEG  +■    1. 

12.  [Get  first   integer  to  be  positive.] 

If  Il(l)*I2(2)    >  0,   go  to  15. 

13.  [Check  node   lengths.]      If  IP1  =    IP2  ,   go  to   Ik  \ 

otherwise,   if  IP1  >   IP2   and  Il(l)    >    0,    go  to 
IT.      Otherwise,   go  to   16. 

Ik. [Check  integer   lengths.]      If  Il(lPl)    >   I2(IP2) 
and  11(1)    >   0,    go  to   17;    otherwise,    go  to  16. 

15. [If  both   <   0,    reverse   sign.]      If  Il(l)    >   0,    go  to  IT. 

16.  [Reverse   signs.]      Il(l)«-   -Il(l),   12  (l)-*-   -I2(l), 
CNEG  «-   -1. 

IT.  [Add  integers.  ]      Il(l+l)  •*■   0.      Il(l)+- 

Il(l)*signum(ll(l))*step(lPl-l)   +   I2(l) 
*signum(l2(l))*step(lP2-l)   +  Il(l).      If 
11(1)    >   BASE,    11(1)  *    11(1)    -   BASE, 
12(1+1)  <-    l.       If  11(1)    <    o,    11(1)  «-    11(1) 
+   BASE,    11(1+1)  «•   -1.      If  I  4  II,   I  «-    1+1, 
go  to  IT. 

18.  [Restore   signs.]      If  Il(l+l)   4  0  ,  I  +■   1+1. 
11(1)  +■   I*CNEG. 
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ALGORITHM  A( A, B,C)  :       C  «-   A+B 

AO.  [Initialize.  ]      BNEG  <-   1. 

Al.  [Find  least    common   denominator.]      Go  to   Ll(A,B,R). 
A2. [A*R,   B*R. ]      Go  to   MO(B,R,T) .      Go  to  M0( A,R,C) . 
A3.  [Convert  to  long  integer.]      Go  to   Cl(C,Il). 

Go  to   Cl(T,I2). 
AU.[Add  integers.]      Go  to  Il(ll,I2). 
A5.  [Convert  back.]      Go  to   Gl(ll,T). 

A6.  [Divide  numerator  by    denominator.]      Go   to  Dl(T,R,C) 
A7. [Restore   B. ]      B(l,l)  +    B(l,l)*BNEG. 


ALGORITHM  S( A, B,C)  :       C  «-   A-B 

SI.  [Initialize.  ]      BNEG  <-    -1,   B(l,l)  «■    -B(l,l). 
S2. [A  +    (-B). ]      Go   to  Al. 

k.      Package   Organization 

Since   the   four  arithmetic   operations  to  be  implemented  all   contain 
two   operands    and  one   result,   setting  each  operator  up   as    a  FORTRAN   subroutine 
with  three   arguments    (essentially  a  three    address   instruction)    is    an 
acceptable    form.      The  user  is   also  provided  with   limited  I/O    capabilities 
involving  other  FORTRAN   subroutines    as    described  in   section   7.      The  two 
conversion   routines    from  section  2    are   additional   routines    called  by  the 
arithmetic  operators.      Since  prime  numbers    are  needed,   there  is   a  routine 
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that  has   some  primes    defined  "by   a  DATA  statement   and  it    can   generate   more 
primes  when   required.      If  a  large  number  of  primes  were   available  on   a 
sequential  access    device,   this  would  be   primarily   an   I/O   file   routine. 

Since  it  is  possible  that  one  of  the  operands  in  an  arithmetic 
operation  may  also  be  the  result,  all  code  is  written  so  that  the  result 
is  not   finalized  until   all  information   is   gotten   from  the  operands. 

An   inspection   of  the  subtraction   and  division   algorithms   reveals 
that  these  need  only  prepare  some   indication  that  the   appropriate   operation 
is  being  done   and  then   use   the  addition   or  multiplication   algorithm, 
respectively.      In   dividing  C  =  A/B  ,   if  IN   is    -1  wherever  C(l,IC)   = 
A(l,IA)   +  B(l,IB)*IN   is  performed,    and  IN  =  +1  when   multiplication  is 
being  performed,   then   the  division  routine  need  only  set   IN  =   -1   and  call 
the  multiplication   routine. 

Similarly,  when  the  subtraction   routine   is   called,  by  setting  an 
indicator  BNEG  =   -1   and  changing  the  sign   of  B(l,l)    for   C  =  A+B ,   subtraction 
is  performed  by   calling  the  addition   routine.      B(l,l)   =  B(l,l)*BNEG   then 
changes   the   operand  back  to   its   original   form  on   exiting  the   addition  routine 

An   alternative  way  to   organize   a  rational   arithmetic  package 
would  be  to   store   two  long  integers    for  each  rational  number  with  tags 
denoting  one   as    a  numerator  and  the  other   as    a  denominator.      The   ordinary 
rules    for  arithmetic  involving  multiple   digit   fractions  would  hold. 
For  example, 

A  *  1  _  k 

B    •   Y        M 
would  be  performed  by 
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L'    =  A*X 
M'    =  B*Y 
R  =   greatest    common   divisor 
(g.c.d.  )    of  L',    M'    found 
by  Euclid's   algorithm   [l] 
L  =  L'/R 
M  =   M'/R 
all  performed  using  routines    for  arithmetic  in  which  both  operands   are 

long  integers . 

A        X        L 
Addition  of  —  +  —  =  —  would  be   accomplished  as    follows. 

L*    =   A*Y  +   B*X 
M'    =   B*Y 
E  =    g.c.d.    of  L'  ,   M' 
L  =   L '  /E 
M  =   M'  /E 
Again,   all  arithmetic  is  performed  on   long  integers. 

If  machine   instructions   are   divided  up  into   sets    according  to 
comparable  execution  times,   then   it   is   possible   to  estimate   the  worst   case 
1   execution  times   of  the   algorithms   needed  with  the  two   different  representa- 
tions   for   rational  numbers.      Denote   the   set   {multiply ,  divide}  by   *,   and 
{ add, subtract }  by  +.      A  long  integer  multiply  with  one  integer  m     digits 
long  and  the   other  m     digits   long  requires   iil  *m     *'s    and,    at  worst,   iil  *nu 
,  +'s.      Using  this,    consider  the   rational  numbers  —with  A  being  n      digits 

D  1 

Y 

long,  B  being  d     digits   long,   and  rr  with   X,Y  having  n    ,    d     digits,   respectively. 
Then   rational   multiplication    and  division   require,    at  worst,   the  number  of 
operations   indicated  in  Table   U.l.      The  value  K   is   the  number  of  divisions 


lU 

required  to   find  the   greatest    common   divisor  R  of  L',   M1  .      Lame    [l]  has 
put  the  upper  bound  on  K  as    ( BASE/2)   while  Knuth    [3]   estimates  K   as 
(12   log  2/iT2)log  N,    for  N   the  minimum  of  L'  ,   M'  .      It  should  he  noted  that 
the  product  of  two  integers  will  have,    at  worst,    as   many  digits    as  the 
sum  of  the  number   of  digits    in   the   operands. 


Multipli  cation 


Reason 


"1  "2 

dl*d2 


K(n1+n2)(d1+d2) 
(n1+n2)*min(n1+n2  ,d1+d2) 
( d]  +  d2 )  *min (n^+rig  A±+^2 ) 


1     2 

Vd2 
K(n1+n2)(d1+d2) 

( n  +n2 ) *min ( n -+n2 , d  + d2 ) 

(dL+d2)*min(n1+n2,d1+d2) 


L'    =   A*X 

M'    =  B*Y 

R  =    g.c.d.    of  L'  ,  M1 

L  =   L'/R 

M  =   M'  /R 


Addition 


n 1*d2+n2*di 

K*max(n1+d2  ,n   +d   ) 
*(d1+d2) 


3/2*max(ni+d2  ,n2+d^)  + 


n1*d2+n2*d2 

Vd2 

K*max(n  +d    ,n  +d   ) 

*(dx+d2) 


L'    =  A*X+B*Y 
M'    =   BY 


min(max(n  +d   ,n  +d   )  ,d  +d   )      min(max(n  +d    ,n  +d   )  ,d  +d   )      L=   L'/E 

*max(  n  +  dg  ,n2+  d   )  *max(  n  +  d2  ,n2+  d1 ) 

min(max(n  +d    ,n  +d   )  ,d  +c^)      min(max(n1+d2  ,n  +d   )  5d1+d2)      M  =   M' /E 


*(d1+d2) 


*(^d2) 


Table   U.l     Worst    case   operation   count  estimate    for  integer 
rational  representation   algorithms. 


15 

In  the  indexed  prime   factorization  representation,   the  only 

multiplications   required  are   due  to   conversion  between   factored  form  and 

long  integers.      Table      k.2   shows  the  worst   case  operation   count  estimate 

for  two   rational   operands   A,  p  nodes   longi   and  B,   q  nodes   long.      The 

value  K1  is  the   average  exponent  of  the  primes,  K        is   the  number  of 

digits   in   the   converted  integer  II,  K       is  that   for  integer  12,    and  K 

"-  3 

is   the  index  of  the   largest   prime   in  the   final  result. 

In  the  sample  program  in  Appendix  I    (using  base    1000   integers), 
average  values    for  the   above  parameters  were   obtained.      These   appear   in 
Table   U.3- 


Multiplication 

* 

+ 

Reason 

p+q 

C(1,I)=A(1,I)+B(1,I) 

Addition 

1/2^*^      +q*K      ) 


l/2K3*max(K21,K22) 


Table    U.2. 


find  g.c.d.    R 

A*R,    B*R 

Convert   A*R-KE1,   B*R-*I2 

add  integers 

convert   to  factored 

form  result. 

Worst   case   operation    count  estimate    for 

factored  form  representation   algorithms. 


p+q 

p+q+2max(p  ,q) 
l/UK1*(p*K21+q*K21) 
3/2   max(K21,K22) 
l/UK3*max(K21,K22) 


16 


PARAMETER 

VALUE 

ni     ■ 

2 

n2 

2 

*1 

2 

^2 

2 

Kl 

2 

K21 

2 

K22 

2 

P 

3 

q. 

3 

K„ 

100 

Table   H.3.      Average   observed  parameter 
values    from  test  program 
A  considerable   savings   in  execution  time   for  the   factored  form 
algorithms   over  the   integer  ration   algorithms    can  be  predicted  on  the  basis 
of  this   test.      Table   k.k  lists  the   comparative    average   values    computed 

using  these  parameters.      Knuth's   estimate   for  K,  based  on  the   average 

6 
integer  being  10   ,  is   approximately  12. 


MACHINE 
OPERATIONS 

INTEGER  RATIO 

FACTORED  FORM 

* 

+ 

* 

+ 

* 

+ 

2  32 
2  32 

2  36 
2^2 

0 

6 

112 
77 

Table   h.k.      Computed  comparisons   using 
test    result  parameters 


B.      RATION— USER'S   MANUAL 

5.      User  Program 

The  RATION  package   currently   available   in  the   PUBLIC  Library 
and  listed  in  Appendix  II  was    designed  to  minimize   -user  responsibility 
for  manipulation   of  rational  variables.      The  user's  primary  responsibilities 
in   the  FORTRAN    calling  program  are  to   dimension  the  variables   properly   and 
to   direct  the  sequence   of  I/O   and  arithmetic   operations   through  FORTRAN 
subroutine   calls.      Several  basic  I/O   facilities    and  all  arithmetic  operations 
are  provided  and  diagnostics    are    available   at  two   levels.      The  list   of  prime 
numbers   is    generated  by  the  package   automatically   although   all  primes    less 
than   1000   are  placed  in  the   array  P   at    compilation  time.      If  subroutine 
RATION   is   called,   a  second  level  of  diagnostics    can    also  be   invoked.      Other- 
wise,   any  routine    can  be    called  with  no   further  preliminary   calls. 

All  rational  variables   used  by  the   routine  must  be   specified  as 
A  is  below: 

INTEGER*2   A(2,K) 
where  K  is   at    least   10.      Rational  arrays    can  be    created  simply  by   adding 
the   array  dimensions    after  the   two   required  dimensions.      For  example, 

B(2,10,5,3) 
is    a  (5,3)    array  of  rational  numbers.      When  these   array  elements    are   used 
as    a  subroutine    calling  argument,   the   (l,l)    element   is  placed  in  the 
argument   list.      For  example,    CALL  ENTER(B(l ,1,2 ,U) )    sets  the    (2,U)    element 
of  B  equal  to   an   integer  read  from  cards.      The   sign   of  any  rational 
variable  is  the  sign   of  its    (l,l)    element;   if  B(l,l)    <   0,   B  is   negative. 

The  program  can  be   initialized  by 
CALL  RATION (N) 
where  N   is   the  value   of  K  if  all   rational  variables    are   dimensioned   (2,K). 
If  during  execution   any  rational  variable  requires    K  larger  than  N,   a 
message   to  that  effect   is  printed.      If  N   is   set  to  0,  no   checking  is   done. 


Any  integers   that   are  not   factorable  in  terms   of  the  available  primes    are 
placed  at  the  end  of  the  prime   array  by  the    conversion   routine  PTOG.      If 
these   integers    (to  the  base  PBASE)    are   more   than  one  word  long,    all 
succeeding  words   of  it    are   tagged  with   a  minus   sign.      When   converting 
back  to   integer  representations,  PI   is   initialized  to  this   unfactorable 
integer  instead  of  to   1.      In  this  way,   rational  numbers   of  reasonable 
complexity   can  be   computed. 

Default    diagnostics    are  printed  when   an   attempt   is   made   to 
divide  by   zero.      The  quotient   is   set  to   zero   and  the  program  continues. 
If  a  conversion  routine   GTOP    (see   section   6)    is    called  to    convert   a 
factored  non-integer  to  integer   form,    a  message  is  printed. 

6.      I/O   and  Utility  Instructions 

Two   input  routines    are  provided  for   creating  factored  form 
integers.      These    can  then  be   divided  to   create   rational  variable  values. 
One   output   routine  writes    a  fixed  format  rational   as    a  signed  two   integer 
quotient.      User   calls    for  the   subroutines    are   listed  below. 

CALL  ENTER(A):      reads    an   integer   from  FORTRAN   FILE   5    (card  reader). 
The   format   of  the   required  cards    follows.      The 
first    card  contains   the   least   integer   greater  than 
or  equal  to  the  number  of   decimal   digits   di viced 
by   four,    right-justified  in   columns   1-5;   -ls  +1  ^n 
columns   6,   7  indicate   the   integer's   sign.      The 
base    10   digits    are   then   divided  in   groups   of  four 
from  the  least   significant   digit  up.      These   groups 
are  placed  starting  with  the  most  significant  group 
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in  each   four  columns   starting  with   column   1. 
For  example,   -177133  has   input 

#00002-l# 

#00177133 #_ 

column     15      8  80 

CALL   CREATE (I,A):       can  be   used  if  the   integer  I    can  "be 

represented  in   an   INTEGER*H   computer  word. 

I    is    a  signed  integer  which  becomes    the  value 

of  A  on   exiting  the   subroutine. 

CALL  WRITE(A):      outputs   A  in   the   form 

+ 1    *   _      X_   _ 
2 

Any   additional   comments   or  spacing  can  be 

provided  by  FORTRAN   I/O  statements   in   the 

user  program. 

If  the  user  feels  he   understands   the   representations    for  rational 
numbers    and  integers    as  presented  in   section   1,  he    can   also   use  two  utility 
routines   that   are    called  by  the  I/O   and  arithmetic   routines.      These  provide 
the   conversions   described  in   section  2.      If  I   is   a  long  integer  variable    and 
A  is    a  factored  form  rational  variable,   then   CALL  PT0G(l,A)    transforms   I   to 
;  A  and  CALL  GT0P(A,l)    transforms   A  to  I .      Thus,    the  program  in   Figure  6.1. 
places   the   factored  form  of  177133  into  A.      The  user   can   save   execution 
time  by  preparing  constants   in  this    and  similar  ways,  but   only  if  he 
understands  what  he   is   doing. 
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INTEGERS   A(2,10) 
INTEGER*  U   1(3) 
DATA  1/3,   7133,  17/ 
CALL   GT0P(I,A) 
Figure   6.1.      Sample   initialization   of 
factored  form  integer 

7.      Arithmetic  Instructions 

Four  arithmetic  routines    are  provided — one   for  each  of  the   four 
operations    closed  on  the  rational  numbers.      RADD   and  RMPY    are  the  principle 
routines  ;   RSUB   and  RDIV  are   supplementary  entries   to  these   using  the   360 
FORTRAN   option   of  multiple  entries    [5].      In   compilers  without  this    feature, 
RSUB   and  RDIV   could  be  very  short   routines   that   call  RADD   and  RMPY, 
respectively,   after  setting  a  special  variable  in   a  COMMON   array.      The 
routines    are    called  as    follows    for  A,   B,    C   factored  form  rational  variables, 

CALL  RMPY(A,B,C):      C   =    A*B 

CALL   RDIV(A,B,C):      C  =   A/B 

CALL  RADD(A,B,C):      C  =   A  +  B 

CALL  RSUB(A,B,C):      C  =  A   -  B 
Any  operators    can  be   redundant. 

CALL   RADD(A,A,A) 
will   double  A,  but  much  more   expensively  than 

CALL  RMPY(TW0,A,A) 
for  TWO  the   factored  form  of  two. 

Logical   comparisons   may  be  made   using  the   fact   that  the  sign  of 
a  rational  variable   is   the  sign   of  its    (l,l)    element.      To   compare   A  with  B, 
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first    compare   the  sign   of  A(l,l)  with  that  of  B(l,l).      If  this    does  not 
clarify  matters  ,  perform 

CALL   RSUB(A,B,TEMP) 
and  compare  the  sign   of  TEMP(l,l).      Obviously,   this   is   quite   time- consuming 
since   RSUB  is   the  most  expensive   routine  to  execute;  however,   comparisons 
are  possible. 

A  sample  program  using  all  arithmetic   routines    and  the  I/O 
routines   other  than  ENTER  is   in  Appendix  I. 
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APPENDIX  I 
Sample  Program 


2k 


IMPLICIT  INTEGEK*2CA-H,  J-0,b)-£) 

IMPLICIT  INTEGEK*4C I ,P) 

COMMON  /COM/  PLIM, PHASE, IPuC  100), PC  1 5000) 

DIMENSION  CC 2, 2 5* 78 >* ON  EC  2, 25),  TMP  (2,  25), 
1  TW0(2,25> 
C 
C************************************************************ 

c 

C  THIS    PKOGKAM    COMPUTES     THE    ELEMENTS     IN    A    LO^EK    TKIANGULAK 

C  12*12    MATKIX.        CCK,M)     IS     KEPKESENTED    AS     CC  K* C K- 1  ) /2  +  M )  . 

C  THIS    PRESENTS    S  i'OKAGE    OF    SPAKSE    ELEMENTS. 

C 

C* **************************************** ****************** 

c 

CALL  KATIGNC25) 
C 

C      C(2,2)=  .5 
C 

CALL  CKEATEC 1,0NE) 

CALL  CKEATEC 2# TWO) 

CALL  KDW  CONE,  TtfG.»C(  1  ,  1  »  3)  ) 
C 

C      C<K»2)=C(K-1*2)  +  1  /  C  2  *  C  K  -  1  >  ) 
C 

DO  10  1= J, 12 

CALL  CKEATEC 1-1 ,  TMP ) 

CALL  KMPY(TWO*TMP*TMP) 

CALL  KDW(ONE*TMH*TMP) 

CALL  KADDCCC \, \,  CI-1 ) *( I -2) /2+2>* TMP*C< 1* 1»I*< 1-1 )/2+2) ) 

CALL  Wril TECCC 1 ,  1 , I*C 1- 1 )/2+2) ) 
10  CONTINUE 

DO  20  1  =  3,  12 
20  CC  1  ,  1 ,  I*C I- 1  )/2+l  )=0 
C 

C        C  C  K  ,  M  )  =  C  C  K  -  1  ,  M  )  +  CM-1)/CM*CK-1))*CCK-1,M-1) 
C  M=3,  .  .  .,K 

C         Ks3#. ..# 12 
C 

DO  30  IK=3,12 

DO  30  IM=3,IK 
CALL  CKEATEC IM*( IK- 1 ), TMP ) 

CALL  CKEATEC IM-1 , ONE) 

CALL  KDW  CONE,  TMP,  TMP  ) 

CALL  KMPYC  CCl,l,CIK-l)*ClK-2)/2+IM-l),  TMP, TMP) 

CALL  KADDCCC l,l,CIK-l)*CIK-2)/2+IM),TMP,CCl,l,IK*CIK-l)/2+IM)) 

CALL  aIKITECCC1,1,1K*C1K-1)/2+IM)) 

IFCPC  l  )  .GE.  1  7 0  00) GO  TO  SI 


*  CCK,  1  ) 


DO  SO  IK=2*12 
KLM=IK*( IK- 1 )/2+ 1 
CALL  CKEATEC  1 ,CC  1 ,  1 ,KLM)  ) 
DO  40  IJ=2,1K 

CCl,l,lK*CIK-l)/2+IJ)  =  CCl,l,IK*CIK-l)/2+IJ)*C-l)**CU+l) 
CALL  KADDCCCl,l,IK*ClK-l)/2+IJ),CCl,l,KLM), 
1  CC  1 ,  1  ,KLM)  ) 
40  CONTINUE 

CALL  /JK1TECCC  1  ,  1  ,KLM  )  ) 

50  CONTINUE 

5 1  S  TOP 
END 


c 

30    CONTINUE 

c 
c 
c 

CCK,  1  )     = 

+     S  UM     C  -  1  )  *  *  1 
1  =  1  ,  .  .  . ,K 
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Q********************************************************************** 

C 
C 
C  *     *     *    KATIGiM     *     *     * 

c 
c 
c********************************************************************** 

c 

C  A    PROGRAM    BY 

C  R.L.     BROWN 

C  MARCH    6,      19  13 

C 

C*  ********************************************************************** 

* 

c 
c 

C    THIS    PACKAGE    READS    AND    CREATES     INTEGERS,     FORMS    RATIONAL    NUMBERS,     ADDS 

C    AND    SUBTRACTS,     MULTIPLIES    AND    DIVIDES    KATIONAL    NUMBERS,     AND    OUTPUTS 

C    RATIONAL    NUMBERS.       ALL    KATIONAL     VARIABLES    ARE    INTEGER*2    AND    SHOULD 

C    BE    DIMENSIONED    CC2,K),     K>10    BUT    USUALLY    LESS     THAT    25.        INTERNAL 

C    REGISTERS    ARE    PLACED     IN    COMMON    BLOCKS    COM,     RT,     P12,     TO     SERVE    SUCH 

C    PURPOSES    AS    STORING     INTEGER*4    PRIME    NUMBERS    AND    HOLDING 

C    INTEGEK    REPRESENTATIONS     IN    FINDING     THE    NUMERATOR     DURING    RATIONAL 

C    ADDITION.     THE    SUBROUTINES    AVAILABLE    ARE: 

C  ENTER(A):     PLACES    AN     INTEGEK     INPUT    FROM    CARDS     INTO    RATIONAL 

C  VARIABLE    A.     SEE    SUBROUTINE    WRITEUP    FOR     CARD    FORMAT. 

C 

C  CREArE(I,A):     PLACES    ANY     INTEGER    OF    LESS     THAT     1  <S     DECIMAL 

C  DIGITS     INTO    RATIONAL     VARIABLE    A. 

C 

C  WRITE(A) : WRITES    THE    CURRENT     VALUE    OF     THE    RATIONAL     VARIABLE    A. 

C 

C  RMPY(A,B,C):     C=A*B  (A,B,C    ARE    RATIONAL     VARIABLES) 

C 

C  KUW(  A,B,C)  :     C  =  A/B 

C 

C  RADD(A,B,C) :     C=A+B 

C 

C  RSUB(A,B,C) :     C=A-B 

C 

C 

C  (GTOP(A,Pl),  PTOG(Pl,A),  GENERP  ARE  CALLABLE  BY  THE  USER  IF  HE  DESIRES 

C      SEE  SUBROUTINE  WRITEUPS  FOR  DETAILS.) 

C 

c*  ******************************************* 


SUBROUTINE    RATION(I) 

C***  ******************************************************  ********* 

c 

C      A    CALL    TO    RATION    WITH     I     SET    TO    SOME    POSITIVE    INTEGER     INITIALISES 

C       A    CHECKING    PROCEDURE    /JH1CH    WAKNS     IF    ANY     FACTORED    FORM    RATIONAL 

C       NUMBER     EXCEEDS     I     IN    LENGTH.       OTHERWISE,     NO     SUCH    CHECKING     IS 

C       DON  E . 

C 

C************************* ********************************  ******** 

IMPLICI T  INTEGER*2(A-H,U-0,U-S) 
IMPLICIT  INTEGEK*4< 1 ,P) 

COMMON  /COM/PLIM,PBASE,PT( 100), P( 10000) 
IF( I ) 1, 1,2 

1  PLIM=2**30-1 
GO  TO  3 

2  PLIM=I 

3  CONTINUE 
RETURN 
END 


27 


SUBROUTINE  GTOPCA,PI) 
IMPLICIT  INTEGER*2CA-H,U-0,U-S) 
IMPLICIT  INTEGER*4C I,P) 

COMMON  /COM/  PLIM,PBASE,PTC 100), PC lOOOO) 
DIMENSION  AC2,  1  ),PI (  1  ) 
C* ****************************************************************** 

c 
c 

C   THIS  ROUTINE  PLACES  INDEXED  PKIME  F ACTOR 1 L AT  ION  FOKM  INTEGER  A  INTO 

C   THE  INTEGER  ARRAY  PI.   PI  IS  INITIALISED  TO  1  AND  MULTIPLIED 

C   SUCCESSIVELY  BY  PRIMES  TO  GET  A.   IE  A  IS  NOT  COMPLETELY  SPECIFIED 

C   BY  THE  AVAILABLE  PRIMES,  THE  UNFACTOkED  PART  IS  TAKEN  FROM  THE 

C   END  OF  THE  PkIME  ARRAY  AMD  INITIALIZES  PI. 

C 

C 

C******************************************************************* 

C  IM     IS    THE    LARGEST    DEFINED    ELEMENT    OF     PI 

C       INITIALISE    PI     TO     1,     IMTIALUE    CONTROL    PARAMETERS. 

IM    =    2 

PIC2)    =     1 

II     =     AC  1,  1  ) 

II     =     IABSCII) 

IF  C I  I  .LE.l  )     GO     TO     12 

PI C  1  )     =     AC  1,  1  )/I I 

IF CAC2, 1 ) .LE.80O0)G0     TO     3 

11=11-1 

1  PI C IM)  =  IABSCPCAC2,  1  )  +  IM-2)  ) 

IF CPCAC2,  1  )  +  IM- 1  )  .GE.O     .OR.     P C A C 2,  1  )  +  1 M- 1  ) . LT . -  9 9999  )  GO     TO    2 
IM=IM+  1 
GO    TO     1 

2  IFC  II  .EO.l  )G0     TO     11 

C       MULTIPLY    MULTIWORD     INTEGER    PI     Br     SINGLE     aIORD     INTEGER    P  C  A  C  2,  I  )  )  *  *  I KK 

3  DO     10    1=2,11 

IF  CAC2,  I )  .GT.PC  1  )  .AND.PC  1  )  .L  T. 8 000) CALL     GENERP 

IF CAC2, I ) .GT  .8000)G0     TO     11 

IKK=AC 1,1) 

IFCIKK.LT.O)GO    TO     13 

DO     10     IK=1,IKK 
IMM     =     IM 

PD=0 

DO     9     IJ=2,1M 

PTC  I J)=PCAC2, I  )  )*PI C IU) 

IF  C IM.EU.2)PD=PTC2)/PBASE 

IFC I  J. EU. 2) GO    TO     5 

PD=PTC IJ-1 ) /PHASE 

PTCIJ-1)=PTCIJ-1 )-PD*PBASE 

PTCIU)     =     PTCIJ)+PD 

PD=PTC I J)/PHASE 
5     IF C IJ .NE.IM.OR.PD. Eu . 0 ) GO     TO    o 

PTC IM)=PTC IM)-PD*PBASE 
A     IMM=IMM+1 

PTC IMM)=PD 

PD=PTC IMM)/PBASE 

IF  CPD.EO.O)     GO    TO    7 

PTC IMM)=PTC IMM )-PD*P BASE 

GO     TO     A 
<>     IF  C  I J  .NE.IM)GO    TO     9 

7  DO    8     ID=2,IMM 

8  PI C ID)=PTC ID) 

9  CONTINUE 
IM= IMM 

10  CONTINUE 

11  PI  C  1  )     =     PI C  1  )*IM 
RETURN 

C       IF    A=<-1,0,+1>    THIS    SECTION    FIXES    PI     TO     THE    PROPER    FORM    FOR    THESE 
C       SPECIAL     VALUES. 

12  PI C  1  )    =    2 
P 1 C  2  )    =     II 

1  K  (  1  1  .N  E  .  0  )     P  1  C  1  )     =     P  I  C  1  )  *  A  C  1  ,  1  ) 
KE  I  URN 
C       ERROR    MESSAGE     IF    A     IS    NOT    AN     INTEGER     CPO/JER    OF    A    PRIME     <    0). 
1    I    WK  ITE(  6,999) 
999    FORMATCSX,  'iriOP    CALLED    /JITH    NON-INTEGER    ARGUMENT') 
I  URN 
F.ND 
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SUBROUTINE  PTOGCPI,A) 

IMPLICIT  INTEGER*2CA-H, J-0,U-Z) 

IMPLICIT  INTEGER*4C I,P) 

COMMON  /COM/  PLIM,PBASE,PTC  1 00 ) > P C  1  0000) 

DIMENSION  PI  (  1  )>A(2j.  1  ) 

DATA  IL/1/ 

C 

c 

C  THIS    KOJTINE    CHANGES        INTEGER    PI      TO     ITS     INDEXED    PRIME    REPRESENTATION 

C  FORM    BY     DIVIDING  BY    SUCCESSIVE    PRIMES    AND    RECORDING    THE    FACTORS. 

C  PI     IS    REDUCED    TO  1     BY     THE    SUCCESSIVE    DIVISIONS. 

C  IE     THE    AVAILAHLE  PRIMES    DON'T    COMPLETELY    FACTOR    PI,     THE    UMFACTOKED 

C  PART     IS    APPENDED  TO     THE    ARRAY     OF    PRIMES. 

C 

C 

c 
c 

C   SAVE  INTEGER  BASE  AND  INITIALIZE  PARAMETERS. 
C 

II  =  IABSCPI CI)) 

IF(PI(2) .LE.  1  .AND.  I  1  .  EU . 2  )  GO  TO  4  5 
3  DO  2  1=1,11 
2  PTCI  )  =  0 

AC  \,   1  )  =  PI ( 1  )/II 
AC l,2)=0 
AC2,2)=0 
IA=2 
IM  =  2 
1  PD=PICII) 
OUT  =  0 
PTC  1  )  =  PI C  1  ) 

c 

C   DIVIDE  BY  PCIM)  AND  SEE  IF  IT'S  A  FACTOR 
C 

DO  10  1=2,11 

IJ  =  I  I  -1  +  2 

PE  =  PD/PCIM) 

PTC  IJ)  =  PE 

IFC I J.EU. IABSCPTC 1 ) ) .AND.PE.EU.O) 
1   PTC  1  )  =  ISIGNC IABSCPTC  1))-1,PTM)) 

PD  =  CPD-PE+PC IM) )*PBASE 
C 

C       SET    OUT    TO     1     IF    PICIM)     IS    NOT    A    FACTOR 
C 

IFCIJ.EU.2.AND.PD.NE.0)     OUT    =     1 

PD    =    PD    +    PI C I  J-  1  ) 
10    CONTINUE 

IF     COUT.EU. 1 )     GO    TO     30 
C 

C       PCIM)     WAS    A    FACTOR; ADD     IT     TO    A 
C 

IFC  IM. EU .AC  2,  I  A)  .OR. AC  2,  1A)  .EU.O)GO     TO     15 

IA=IA+1 

IFC IA.GT  .PLIM)GO     TO     4  1 
AC  1,  IA)  =  0 
15    AC  1,  IA)=AC  1,  IA)+ 1 

AC  2,  1A)  =  IM 

I  I     =     IABSCPTC  1  ) ) 

DO    20    IJ=  1  , I  I 
20    PI C IJ)     =    PTC IJ) 

IF C I  I  .EU. 2 .AND. PI C2)  . EU .  1  )     GO     TO     40 

GO    TO     1 
30     IF  C  IM+  1  .GT.PC  1  )  .AND. PC  1  )  .LI".«000)     CALL     GENERP 

I M  =  I M  +  1 

IF C IM.LE.HOOO)GO    TO     1 

IFCAC2,IA).NE.0)IA=1A+1 

AC2, IA)=HOOO+IL 

AC  1,  IA)=  1 

DO     .)',     1  =  2*11 
35    PC8000+I-2+lL)=PICI)*iSlONCl,5-2*I) 
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IL=1L+II-1 

PCIL+8000)=-9999999 

GO    TO     40 

41 

WRITEC6*999) 

999 

F0RMATC5X,  'FACTORED 

40 

AC  1 >  1  )=IA*AC  1,  1  ) 
AC2, 1 )=AC2, IA) 

42 

RETURN 

c 

c 

IF 

PI     =     <-l,0,+l>,     SET 

c 

FORM     VARIABLE    NOT    LARGE    ENOUGH*) 


A    TO     ITS       SPECIAL    FORM. 


4  5 


AC  1,  1  )     = 
AC2,1)     = 
RETURN 
END 


PI  C2)*PI C  1 )/I I 
1 


THIS    ROUTINE    GENERATES     THE    NEAT    PRIME    IN    THE    ARKAY    MPCI)     JSING 

KNUTH'S    ALGORITHM    P,     SECTION    1.3-2 

MPCM+1)     CONTAINS    A    POSSIBLE    PRIME    NUMBER. 

MPU     IS    THE    UUOTIENT    OF    DIVIDING    MP(M+1)     BY    A    KNOWN    PRIME; 

PR    IS    THE    REMAINDER. 


SUBROUTINE    GENERP 

IMPLICIT     INTEGER*2CA-L, N-^) 

REAL*8    PR 

COMMON    /COM/    MLIM, MBASE, MTC 100), MPC 10000) 
C 
C*********************************************************  ********* 

c 
c 
c 
c 
c 
c 
c 
c 
c 

c*  *  *  *  *  *  ************************************************************* 

c 

M    =    MP (  1  ) 

MPCM+1  )    =    MPCM)    +    2 
1     DO     10    1*2* M 

MPU    =    MPCM+  1  )/MPC I  ) 

IF CMODCMPCM+ 1  ),MPC I  )  )  .  EU  .0)G0    TO     12 

IFCMPQ.LE.MPC I ) )    GO    TO     11 

10  CONTINUE 

11  MPC  1 )  =  M+  1 
RETURN 

12  MPCM  +  1  )=MPCM+ 1  )+2 
GO  TO  1 

END 


BLOCK  DA 

COMMON  / 

INTEGER* 

INTEGER* 

+53*59,61 

+137, 139, 

+  223,227, 

+307, 31 1, 

+389, 397, 

+479, 487, 

+  577,  587, 

+  659,661, 

+757, 761, 

+857,859, 

+953,967, 

END 


TA 

COM/  ML  I 
4  MBASE/ 
4  MP/ 169 
,67,71,7 
149, 151, 
229,233, 
31 3, 317, 
401, 409, 
49 1 , 499, 
593*599* 
673,677, 
769, 773, 
863,877, 
971  ,977, 


M, MBASE, 
10000/, M 
,2, 3, b,  / 
3,79,83, 
157,  16  3, 
239,241, 
331, 337, 
419,421, 
503, 509, 
601,607, 
68  3,6  9  1, 
787, 797, 
881  ,883, 
983,991, 


MTC  1 
LIM/ 
,11, 
89,9 
167, 
251, 
347, 
431, 
52  1, 
613, 
701, 
809, 
887, 
997, 


00), MPC 10000 
9999999/ 
1  3,  1  7,  19,23, 
7,101, 103, 10 
173, 1 79, 181, 
257,263,269, 
349, 353, 359, 
433, 439, 4$3, 
523, 541, 547, 
6  17, S  19,6  31, 
709, 7 19, 727, 
8  1 1,82  1,823, 
907,91 1,919, 
9831*0/ 


) 

29,3 
7,  10 
19  1, 
271, 
36  7, 
449, 
557, 
641, 
733, 
827, 
929, 


1,  37 
9,  1  1 
193, 
277, 
37  3, 
457, 
56  3, 
643, 
739, 
829, 
937, 


,41,  43,  47, 
3^  127,  131, 
197,  199,21  1 j 
281,283,293; 
379, 383, 
46  1,  46  3; 
5  69, 57  1, 
647,6  53, 
743,  751, 
839,853, 
941,947, 


467i 
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SUBROUTINE    KMPYCA,B,C) 
IMPLICIT     INTEGEK*2CA-H, J-0,U-Z) 
IMPLICIT     INTEGER*4C I,P) 
INTEGER* 4    MINCMAXO 

COMMON    /COM/    PLIM,PBASE,PTC  100), PC  10000) 
DIMENSION    C(2* 1)*AC2* 1)»B(2* 1)»D«2»50) 

C       THIS    ROUTINE    MULTIPLIES    OK     DIVIDES    BY    ADDING    OR    SUBTRACTING 

C       THE    EXPONENTS    OF    THE    PRIME    NUMBERS     IN    THE    FACTORIZATION    OF    THE 

C       OPERANDS.        IF    A    PRIME    IS     DEFINED    FOR    A    BUT    NOT    FOR    B,     ITS     EXPONENT 

C       IN    B     IS    0.        ENTRIES    ARE: 

C  RMPY:     A*B 

C  KDIV:     A/B 

C*  ************************************  **************************** 

C    INITIALIZE    PARAMETERS 

C 

IN    =     1 
1     PA=AC1,1) 

PB    =    BC 1, 1 ) 

PA=IABS(PA) 

PB=IABSCPB) 

IF (PA*PB.LE.l )     GO    TO     25 

ID=2 

IA  =  AC2,  1  ) 

IB=BC2» 1 ) 

II=MAX0C IA, IB) 

IA=MIN0CPA,2) 

IB=MIN0CPB,2) 

DC2, 1 )=1 

DO     10    1=2,11 

IF C AC  2, I  A) *1  .EO. I )G0     TO     3 

GO    TO     4 

3  IF CAC2, IA)  .EU.BC2, IB) ) GO     TO    ^ 
DC 1,  ID)=AC 1, IA) 

DC2, ID)=AC2, IA) 
IA=MINOC IA+l.PA) 
GO    TO    9 

4  IFCBC2, IB)*1  .NE.I ) GO     TO     10 
DC  1, ID)  =  BC  1, IB)*IN 

DC2, ID)=BC2,IB) 

IB=MIN0CIB+1,PB) 

GO    TO    9 
6     DC 1, ID)=AC 1, IA)+IN*BC 1, IB) 

DC2, ID)=BC2, IB) 

IA=MINOC IA+ 1,PA) 

IB=MINOC IB+ 1,PB) 
9     IFC I .  EU. I  I )G0    TO     10 

IF  CDC  1, ID)  .NE.O)     I D= I D+ 1 

IF  CID.GT.PLIM)GO    TO     13 
10    CONTINUE 
C       D    IS    TEMPORARY     STORAGE    SINCE    A    OR    B    MAY     EUUAL    C    IN    THE    CALLING 
C       ROUTINE. 

IF  CDC  1,  ID)  .EU.O)ID=IU-  1 

DC  1,  1>  =  ISIGNCIDjA<  1,  1 )*BC  1,  1 )* 1 ) 

DC2, 1 )=DC2, ID) 

GO    TO     \A 

13  WRITEC6,999) 

999    FORMATC 5X, 'FACTORED    FORM     VARIABLE    NOT    LARGE    ENOUGH') 

14  DO     15    1=  1 ,  I D 
CC 1, I )=DC 1,1) 

1 5  CC2, I )=DC2, I ) 
RETURN 

1,0,  +  1>,     COMPUTE    THE    SPECIAL     FORM    OF    C 
PRINT    AN    ERROR    MESSAGE. 


C 
C 

c 
c 

IF     BOTH    A 
IF    B    IS    A 

AND    B    ARE    IN    <- 
ZERO     DIVISOR,     P 

25    CC 1, 1 ) 

=    AC  1,  1  )*BC  1,1) 

CC2,  1  )=1 

IFCPB.EO.U     .AND.     IN.EU.-I)     WK I TE Co , 1 000) 
1000    F0KMATC5X, 'DI VISION    BY     ZEROJ     QUOTIENT    SET    TO    ZERO') 
RETURN 

ENTRY    R[)IVCA,B,C) 
IN    =     -1 
GO    TO     1 
END 
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SUBROUTINE    KADD(A,B*C) 

IMPLICIT    INTEGER*2(A-H, J-G,u-^) 

IMPLICIT    INTEGER*4( I,P) 

INTEGER*4    MINO,MAXO 

COMMON    /COM/    PLIM,PBASE,PTC 100), P( lOOOO) 

COMMON     /KT/R(2,50),T(2, 50) 

DIMENSION    AC2,  1),B(2,  1  ),CC2,  1 ) 

COMMON    /P12/    P1C50)*P2(50) 
C 
C 

c 
c 

C  THIS    ROUTINE    ADDS    OK    SUBTRACTS    A    AND    B,     PLACING    THE    RESULTS     INTO     C 

C  ENTRIES    ARE: 

C  KADDCA,B,C)  2     A+B=C 

C  RSU8(A,B,C):     A-B=C 

C 

C 

Q+  if  if  if  if  ^c  ^c  if  :+:  if  if  if  if  if  ■%.  if  if  if  if  if  if  if  -if.  if  if  if  if  if  if  if  if  if  if  if  if  ^c  +%^c*:  +  ^t  +  ++5f:^t^;^  +  +  +  ^:+:^;+3f:+^c%^:^c^:+:^a|c  + 

c 

C       PF(I,P)     IS    USED     IN    ADDING    TWO     INTEGEKS    OF    DIFFERENT    LENGTH 

C       BY    MERGING    THE    SHORTER    WHEN     NECESSARY. 

C 

PFC  IZ,PC)=MAXO(PC-I^+  1jO)/MAXO(PC-I^+  l>l) 

BNEG    =     l 

GO    TO     l 

ENTRY    RSUBCA,B,C) 

BNEG    =    -l 

BC  l,  l )     -    -B(  l,  l  ) 

1  PA  =  AC  l,  l  ) 
PB  =  B( l , l ) 
PA=IABS(PA) 
PB=IABS(PB) 

IF( < AC  2, PA- l  )  . GT.8000.AN0.PA.NE.2)  .OR . ( BC 2, PB-  l  ) .GT.8000 
l      .AND.PB.NE.2) )GG    TO     llU 
CNEG    =     l 
DO    2    I = l , l 00 
KC  l  ,  I  )  =  0 

2  PT(I )    =    0 
RC  l,2)=0 
R(2, l)=l 
IR  =  2 

II=MAX0(A(2, l )*1#B<2» 1 )*1) 

IA=MIN0(2,PA) 

IB=MIN0(2,PB> 

IFUI  .LE.  1  )G0    TO     13 
C 

C      FIND    THE    LEAST    COMMON    DENOMINATOR    BY     STOKING    THE    MOST    NEGATIVE 
C       OF    ALL     NEGATIVE       EXPONENTS    OF    PRIMES. 
C 

DO     10    1=2,11 

IF(A(2,IA) .EG. I )G0    TO    3 

GO    TO     A 

3  IFCAC2, IA) .EU.BC2,  IB))GG    TO    6 
IF (AC  1, IA) .GE.O)GG    TO     5 

R(  1,  IR)=A(  1* IA) 
H<2, IR)=A(2, IA) 
IR  =  MINO( IR+1, I  I  ) 

5  IA=MINO( IA+1,PA) 
GO    TO     10 

A    1F(B(2, IB)  .NE.  I  ) GO    TO     10 

IF CB( 1, IB) .GE.O)    GO    TO    7 

R(  It IK)=B<  1, IB) 

K<2, IR)=BC2, IB) 

IR=MINO( IR+ 1*11) 
7     IB=MINO< IB+ 1,PB) 

GO    TO     1U 

6  R(  1, IK)=MIN0(  1*A<  1, I  A),  1*B(  1, IB) ) 
IF(R(  1,  IR)  .QE.OGO    TO    9 
R<2»IH)=A(2#IA) 

IR  =  MINO( IR+ 1, I  I ) 
9     IA=MINO( IA+ 1,PA) 
IB=MINO( 18+1, PB) 
10    CONTINUE 
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C 


IF(KC 1 

,  i  r ) . 

GE. 

0)IR= 

IK- 

■  1 

KC  1 

,  1  ) 

=  IR 

Hi  2 

,  1  ) 

=  K(2, 

IK) 

GO 

ro 

1  4 

K<  1 

*  i  > 

=  1 

R(2 

,  i  ) 

=  1 

13 

C 

C   PUT  PARTIAL  NUMERATORS  INTO  T  AND  C  BY  A*K,  B*R,  WHERE  R  IS  THE 

C   LEAST  COMMON  DENOMINATOR.   CONVERT  TO  INTEGERS  P1,P2. 

C 

14    CALL    RDIV(B,K,T) 
CALL    RDIV(A,K,C) 
CALL    GT0P(C*P1) 
CALL    GTGP(T*P2) 
C 

C       ADD    P2     rO    PI     AFTER    ADJUSTING    SIGNS     TO       GET    A    POSITIVE    RESULT. 
C 

IP1    =    IABSCP1 ( 1 ) ) 
IP2    =    IABS(P2( 1 ) ) 
IF(P1 ( 1 )*P2( 1 )  .GT.O)     GO     TO    7  0 
IFUP1-IP2)     30,20*40 
20    IFCPK  IP1  )  .GT.P2(  IP2)  )     GO    TO     40 
30    IFCP1 ( 1 )  .LT.O)     GO    TO    60 
35    PI ( 1 )     =    -P 1 ( 1 ) 
P2( 1 )     =    -P2(  1  ) 
CNEG    =     -  1 
GO    TO    8  0 
40    IF(P2( 1 ) .LT.O)     GO    TO    80 

GO    TO    3  5 
70    IF (PI  (  1  )  .LT.O)     GO    TO    35 
80    II     =    MAXOC IP1, IP2) 
PT(2)=0 
DO    90    1=2, I  I 
PT(I+1 )=0 

PTCI)    =    P 1 ( 1 )/IPl *P1 ( I )*PF( I, IP1 )+PF( 1, IP2)*P2( 1 )/IP2*P2( I ) 
1     +PTCI) 
IF    (PTC  I ) .GE.O)    GO    TO    8  5 
PT(I)    =    PT(I)    +    PBASE 
P  T  (  I  +  1  )    =    PT(I+1)       -     1 
GO    TO    9  0 
85    IF(PT( I ).LT. PBASE)     GO    TO    90 
PT(I)    =       PT(I)    -    PBASE 
PT( 1+ 1 )    =    PT( 1+1 )     +     1 
90    CONTINUE 

IF(PT( I  1  +  1  ) .NE.O)     11    =     11     +     1 
DO     100    1=2,11 
100    PI  ( I )    =    PT( I ) 

Pl(l)    =    II     *    CNEG 
C 

C       CONVERT    FINAL    NUMERATOR     TO     INDEXED    PKIME    FACTORISATION    FORM, 
C       DIVIDE    BY     K,     AND    PLACE    ANSWER     IN    C. 


CALL    PT0G(P1,T) 
B(  1,  1  )     =     B(  1,  1  )*BNEG 
CALL    RMPY(T,K,C) 
RETURN 
1  10    a)RITE(6,  1000) 
1000    F0RMAT(5X, 'ADDITION    OPERAND    HAS    MORE    THAN    ONE    UNFACTORED 
1, 'COMPONENT;     RADD     JILL    YIELD     INCORRECT    ANSWER.') 
RETURN 
END 
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C* 

C 

C 

C 

C 

C 

C 

C* 


SUBROUTINE  WKITECA) 

IMPLICIT  INTEGER*2CA-H,J-0,U-/) 

IMPLICIT    INTE6ER*4C  I*P) 

INTEGEK*4    MAXO>MINO 

COMMON    /COM/    PLIM,PBASE.»PTC 100), PC lOOOO) 

COMMON    /RT/RC2,50),TC2, 50) 

DATA     IDASH/' '/ 

COMMON     /P12/P1 C50)>P2C50) 

DIMENSION    AC  2,  1 ), IFMC30), I  FORM (4  )  >  IF  1  (  1  2 ) > I F2(  1 2 ) , I F 3(S  ) 
EQUIVALENCE    C I FM C 1  )  ,  IF  1 ( 1 ) ) , ( I FM (  1 3 ) , IF2C 1 ) ) , C I FM ( 2 5 ) , I F3C  1 ) ) 
DATA    IF1/' C  1H+,23X,27I4)  C  I  H  +  ,  39X  ,  221  4, /7X,  301  4  )     '/ 

DATA    IF2/' C 1H+,55X, 19I4,/7X, 3014)     (  1 ri+ ,1  IX ,  1 51 4, /7X ,  301 4 )     V 
DATA    IF  3/' ( 1H  +  ,87X,  101 4,/7X, 3014)      '/ 


THIS    ROUTINE    WRITES    OUT    TWO     INTEGERS 
TO    REPRESENT    A    BASE     10    RATIONAL     NUMBER 


II     =    AC  1,  1  ) 

II     =     IABSCII) 

IF  CII  .LE.  1  )     GO     TO    26 

IA  =  2 

IR  =  2 

IT  =  2 


PLACE    THE    NUMERATOR    (POSITIVE    EXPONENTS) 
(NEGATIVE    EXPONENTS)     IN    K. 

RC2,  1  )=1 
TC2,  1  )=1 

CONVERT    BOTH    TO     INTEGERS    AND    WRITE. 


2     IFCAC  }j  IA)  )4,4,6 
4    R(2, IR)=A(2, IA) 

RC 1, IR)=-A( 1  ,  IA) 

IK=IR+ 1 

GO    TO     7 

6  TC2, IT)=AC2> IA) 
TC  \, IT)=AC  1, IA) 

IT=IT+1 

7  IA=IA+1 

IF  (  IA.LE.I  I  )G0    TO    2 

RC  1,  1  )  =  IK-1 

RC2>  1  )=RC2, IR-  1  ) 

TC  1,  1  )    =    IT-1 

TC2,  1 )=TC2, IT-1  ) 

25  II    =    AC  1,  1 )/I I 
IK=IR-1 

I T= IT-1 
GO    TO    29 

26  TC  1,  1 )  =  AC  1,1) 
TC2>  1  )  =  1 

R(  1,  1  )=1 
RC2, 1  )=  1 

29  CALL    GTGPCTjPI ) 
IS  =  0 

30  IFCTC2, IT- 1 ) .LE.HOOO.OR.IT. EU . 2  )  GO     TO    33 
IT=IT-1 

IS=IS+  1 

ISS=-1 

PTCIS)=99999 

31  IS=IS+1 
ISS=ISS+1 

PTCIS)=IABSCPCTC2, IT)  +  ISS)  ) 
IFCPCTC2, IT)  +  ISS+1  )  .L  T.O.AND. 

1     -9999999)G0     TO     31 
GO    TO    30 
33    IS=IS+1 


IN    TJ     DENOMINATOR 


PCTC2,  ID  +  ISS+  1  )  .GT. 


3k 
PTCIS)=99999 
I  1  =  IABS(PJ C 1 ))-l 
IFR  =    MGD<  I  l,20)/4 
DO    9  9     10=1,6 
99     IFGKM<IG)=IFMC IG+6*IFK) 

•JK I  TEC  6,  1001  )     (Pl(Il+2-I)>I  =  l.Il) 

IFCIS.GT.l )WKlTEC6>  I  FORM )C PTC  IS+2- I ) > I =2* I S) 

I'S=0 

CALL  GTGPCk,P2) 

34  IF  CRC2,  Ir<-  1  )  .LE.8000.GK  .IT.  Eu  .2)GG    TO    37 
IK=IK-1 

IS=IS+1 

I S  S~- 1 

PTOS)=99999 

35  IS=I$+1 
ISS=ISS+1 
PTCIS)=PCRC2, IR)+ISS) 

IFCPCRC2,IR)+ISS+1)  .LT.G.AND.PCKC2, IiO  +  ISS+1  )  .GT. 
1     -9999999) GG    TG    35 

GG    TG    34 
37     I S=  1S+  1 

PTCIS)=99999 

I2=IABSCP2C  1  )  )-l 

IJ=MINOCMAXOC I  1 , 12)* 20) 

WRITE (A*  1002)     II, C I  DASH, IK=  1, IJ) 

IFR  =    M0DCI2,20>/4 

DO    98     10=1,6 
98     IFOKMC IG)=IFMC IG+6*IFR) 

WK I TEC6>1001)CP2CI2+2-I), 1=1,12) 

IFCIS.GT.I)  WKITEC6, 1F0KM)  C PT C I  5  +  2- I ) , I  =  2, I  5 ) 

RETURN 

1001  FOKMATC 5( 6X,  SC  IX, 201  4/) ) ) 

1002  FGKMAT(2X, 12, *    *     ',20A4) 
END 


SUBROUTINE    CREATECI,A) 

IMPLICIT    INTEGER*2CA-ri, J-0,U-<O 

IMPLICIT    INTEGER*4C I,P) 

COMMON    /P12/PI C50),P2C SO) 

DIMENSION    AC  1,  1  ) 

COMMON    /COM/    PLIM, PHASE, PTC 10100) 
C 
C*  *************************************************************** 

c 
c 

C      THIS    ROUTINE    CREATES    THE    INDEXED    PRIME    REPRESENTATION    FORM    OF 
C      ANY    INTEGER    I    REPRESENTABLE    Hi    THE    MACHINE.       PI     IS    REDUCED    TO     1 
C 
C 

C*  ********************  ************************************  ****** 

c 

IM    =    2 

II     =     IABS(I) 

IF  <  1 1  .LE.  1  )    GO    TO    20 

PI  (  1  )    =     I/I  I 
1     PI ( IM)     =    MODCI I,PBASE) 

II     =     CII-PK  IM)  )/HHASE 

IF ( II .EO.O)GO    TO     10 

IM    =     IM+1 

GO    TO     1 
10    PI (  1 )    =    PI  (  1  )*IM 

CALL    PTGG(PI,A) 

RETURN 
20    PI ( 1 )    =    2 

PIC2)    =    I  I 

IF(II.NE.O)     Pl(l)     =    P  I  C  1)  *  1 

AC  1,  1  )     =     I 

A  <  2  ,  1  )  =  1 

RETURN 

END 


35 


SUBROUTINE    ENTER(A) 

IMPLICIT    INTEGER*2CA-H, J-0,U-Z) 

IMPLICIT    INTEGER*4C I,P) 

COMMON    /COM/    PLIM,PBASE,PTC 100), PR( 10000) 

COMMON    /P12/    P(50),P2(50) 
C 
C***************************************************************** 

c 
c 

C  THIS  ROUTINE  READS  LONG  INTEGERS  FROM  CARDS  AND  CONVERTS  THEM 

C  TO  A.   PI  IS  SET  TO  1.   THE  FORMAT  IS 

C  FIRST  CARD: 

C  IN=GREATEST    INTEGER    PART    OF    (  #GF    DIGITS)/4    IN    COL.     1-5; 

C  SIGN    =<-l,+ 1 >    IN    COLUMNS    6,7. 

C  FOLLOWING    CARDS: 

C  DIGITS    OF    THE    INTEGER    IN    GROUPS    OF    FOUR    DIGITS,    RIGHT    JUSTIFIED 

C  WITHIN    GROUPS,     GROUPS    LEFT    JUSTIFIED. 

C  EXAMPLE:       -177133 

C  CARD    ONE    #00002- 1# 

C  CARD    TWO    #001771 330000.. .0000# 

C 

C 

C* **************************************************************** 

READC5, 1000)  IN,ISIGN 

READC5, 1001 )     (PC IN+2-I ), 1= 1, IN) 

P(  1  )    =     CIN+1  )*ISIGN 

CALL    PTOG<P,A) 

1000  FORMATC 15,  12) 

1001  FORMAK  5(201  4)  ) 
RETURN 

END 
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